Simulation of melting of two dimensional Lennard-Jones solids 
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We study the nature of melting of a two dimensional (2D) Lennard-Jones solid using large scale 
Monte Carlo simulation. We use systems of up to 102,400 particles to capture the decay of the 
correlation functions associated with translational order (TO) as well as the bond-orientational 
(BO) order. We study the role of dislocations and disclinations and their distribution functions. 
We computed the temperature dependence of the second moment of the TO order parameter ($g) 
as well as of the order parameter ty^ associated with BO order. Applying finite-size scaling of 
these second moments we determined the anomalous dimension critical exponents n(T) and r]e(T) 
associated with power-law decay of the and ^6 correlation functions. We also computed the 
temperature dependent distribution of the order parameters ^g and on the complex plane which 
support a two stage melting with a hexatic phase as an intermediate phase. From the correlation 
functions of ^g and we extracted the corresponding temperature dependent correlation lengths 
£(T) and ^6{T). The analysis of our results leads to a consistent picture strongly supporting a 
two stage melting scenario as predicted by the Kosterlitz, Thouless, Halperin, Nelson, and Young 
(KTHNY) theory where melting occurs via two continuous phase transitions, first from solid to a 
hexatic fluid at temperature T m , and then from the hexatic fluid to an isotropic fluid at a critical 
temperature Ti. We find that £(T) and £e(T) have a distinctly different temperature dependence 
each diverging at different temperature and that their finite size scaling properties are consistent 
with the KTHNY theory. We also used the temperature dependence of rj and 776 and their theoretical 
bounds to provide estimates for the critical temperatures T m and Ti, which can also be estimated 
using the Binder ratio. Our results are within error bars the same as those extracted from the 
divergence of the correlation lengths. 
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PACS numbers: 64.60.De,67.70.D-,61.72.Bb 



I. INTRODUCTION 

The most widely considered theory of 2D melting is 
the so-called KTHNY theory of Kosterlitz and Thouless 1 , 
Halperin and Nelson 2 ^, and Youngi, which predicts that 
melting in two dimensions occurs via two continuous 
phase transitions, first from solid to hexatic fluid, and 
then from hexatic fluid to isotropic fluid. This theory 
begins from the fact that true translational order can- 
not exist at any non-zero temperature in 2D because of 
the infrared divergence caused by the zero point motion 
of long-wave-length density fluctuations. According to 
the KTHNY theory, another form of true long-range or- 
der exists below some non-zero temperature T m where 
only the directions of the nearest-neighbor bonds order. 
This long-range bond order disappears above T m because 
of dislocation unbinding which leads to an intermediate 
phase, the hexatic phase, characterized by topological or- 
der, where while dislocations are unbound, disclinations 
with opposite topological charge remain bound. These 
disclinations become unbound at a higher temperature Tj 
where the system becomes an isotropic disordered fluid. 

Simulation of melting in classical two-dimensional (2D) 
systems has been tackled by means of a variety of com- 
putational studies^ for several decades without reaching 
a definite conclusion regarding its nature. In particular 
for hard disks in 2D, a large number of computer simula- 



tion studies have been applied to understand 2D melting, 
since this is the toy model on which the Metropolis Monte 
Carlo method itself was first introduced^ and soon after- 
ward, the 2D melting of hard disks was studied 7 . One 
of the reasons for the difficulty to reach an unequivocal 
conclusion is that in 2D a conventional solid with true 
translational order cannot exist, and, instead the corre- 
lations decay very slowly over long distance. This re- 
quires large size systems where the relaxation time scales 
become very long for these types of phenomena. In par- 
ticular for hard disk systems, when using a local updating 
algorithm or even molecular dynamics, particles remain 
stuck in their local "cage" for large computational time 
scales, precisely because of the hard disk constraint. 

One might think that Monte Carlo simulation of soft- 
core potentials, such as the Lennard-Jones system in 2D, 
might not be plagued by the same level of computa- 
tional severity as the hard-disk systems, because of the 
softening of the hard-core constraint. As a matter of 
fact there are a number of studies of the Lennard-Jones 
solid 8 by computer simulation where also a general con- 
sensus about the nature of melting has not been estab- 
lished. Some studies have favored a first-order transition 
from solid to liquid^—, as predicted by the grain bound- 
ary melting suggestion 1 —, while other studies^ - — have 
leaned toward the KTHNY theory. The most thorough 
of these studies, however, are at least one decade old and 
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because of the fact that the computational resource con- 
straints of today are significantly better, a more thorough 
study should be possible. 

In the present paper, we study the nature of melt- 
ing of a two dimensional (2D) Lennard- Jones solid using 
large scale Monte Carlo simulation. We use systems of 
up to 102,400 particles to capture the decay of the corre- 
lation functions associated with translational as well as 
the bond-orientational order. We find that to carry out 
thorough investigations beyond these sizes, calculations 
using the Metropolis local update become impractical us- 
ing today's high performance computing because of the 
long relaxation time scales. Further technical details of 
our simulation arc described in the next section, and the 
remainder of the paper is organized as follows. In Sec. lIIII 
we discuss the role of defects in the KTHNY theory of 
melting and present the results of a geometric defect anal- 
ysis. In Sec. IIVI we show the temperature dependence 
of both order parameters, an d as well as their 
second moments, and ^>\. The system-size depen- 
dence of and can be used to determine the criti- 
cal exponents r\ and rje, as shown in Sec. [V] In the same 
section, the KTHNY values of the critical exponents at 
melting, rj(T m ) and ?yg(Ti), are used to estimate the tran- 
sition temperatures T m and X^. Next, in Sec. IVI[ we 
present our results on the correlation function associated 
with bond orientational order above Tj and determine 
the temperature dependent correlation length £e(T). In 
the same section, we demonstrate finite-size scaling of 
\&g. A similar presentation is given in Sec. IVIII for the 
pair distribution function and the correlation length of 
translational order, £(T). In addition, we present our 
findings for the scaling behavior of the second moment 
of in this same section. Sec lVIUl presents an analysis 
of the melting transition using Binder's cumulant ratioi^ 
for each order parameter, and also includes a discussion 
of finite-size scaling theory in the presence of multiple 
correlation lengths. Finally, in Sec. IIX1 we give a brief 
summary of our main findings and conclusions. 

II. SIMULATION DETAILS 

In the Lennard- Jones potential, for two particles sep- 
arated by a distance r, 

™ -KG)" -(?)")■ (1) 

an attractive inverse sixth power tail is combined with a 
repulsive inverse twelfth power hard core, such that there 
are only two parameters: e, the potential well depth, and 
<7, the hard-core diameter. However, our results can be 
inferred for any particular value of these parameters (for 
a specific real system), since in our calculations, distance 
is measured in units of a and temperature in units of 
e/k B . 

In our calculations we have truncated the Lennard- 
Jones potential at a distance of 3a, and shifted the value 



of the potential within this cutoff distance by a con- 
stant so that the resulting potential approaches zero at 
3er (matching the values beyond 3cr). This truncation 
is justified because the Lennard- Jones potential is al- 
ready quite small (-0.005e) at this distance, and is not 
expected to significantly affect the accuracy of our simu- 
lations. Additionally, by using a cutoff distance, we are 
able to use a cell list structure in our algorithms so that 
our computations scale as 0(N) instead of the 0(N 2 ) 
scaling without a cutoff distance (N is the number of 
particles in our simulation cell). 

We have collected data for systems of 1600, 6400, and 
25600 particles over a wide temperature range at a den- 
sity of 0.873. Additionally, we have simulated a system 
of 102400 particles for three temperature values at the 
same density in order to verify our results for the smaller 
system sizes. To accommodate the expected low temper- 
ature triangular solid phase, a periodic simulation cell 
of proportion 2:\/2 is used. We have performed our cal- 
culations on the Florida State University shared High- 
Performance Computing facility, which contains several 
thousand compute nodes. The processors on these nodes 
range in speed from 2.3 GHz to 2.8 GHz, and it takes 
about 34 hours to perform 1,000,000 Monte Carlo sweeps 
for N — 25600 particles, including calculating observ- 
ables every 100 Monte Carlo sweeps (MCS). We have 
found that, except for the N = 102,400 particle system, 
a million MCS are sufficient to reach equilibrium, even 
near the critical points. The data presented here is ob- 
tained over one million MCS, after a period of one (for 
N = 102,400) or two (for N = 1600,6400, and 25600) 
million MCS of equilibration. 

To take advantage of our computational resources, we 
utilized a trivially parallel Monte Carlo implementation 
of 100 threads, each with a unique random number seed 
and initial configuration. Simulations begin from an ini- 
tial near-ordered configuration (particles are placed in a 
triangular lattice, with 5% lattice spacing random fluctu- 
ations). Statistics for thermodynamic variables are col- 
lected by generating averages on each of the 100 parallel 
threads, then using the central limit theorem we obtain 
the total average, as we have 100 independent means. 

Although in our preliminary studies we have computed 
thermodynamic quantities for a range of densities and 
temperatures, the effects of critical slowing down near 
the melting transition and our desire to study the largest 
possible systems have led us to focus on a single den- 
sity, 0.873 (all densities are in units of particles per cr -2 ). 
This density was chosen for several reasons. This is a 
density that could be readily compared to prior numer- 
ical simulations of Lennard- Jones melting 15 . Also, we 
wanted a density that is relatively low, but large enough 
to avoid the solid-vapor coexistence phase at low tem- 
peratures. Strictly speaking, there is a solid phase in the 
zero temperature limit only at densities of 0.9165 (the 
density at which the spacing of the triangular lattice is 
the same as the position of the Lennard- Jones poten- 
tial minimum) and above. Below this density there is a 
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solid-vapor coexistence phase. However, the triple point 
density is roughly 0.82, so at higher densities the system 
will in general become solid before the onset of melting 
occurs^. 



III. ROLE OF DEFECTS 

A. Defect types 

In two dimensions, the densest packing of particles of 
uniform size is achieved in a triangular lattice. In such a 
configuration, each particle has exactly six nearest neigh- 
bors. Thermal fluctuations will lead to distortions in the 
lattice, or even destroy it completely. To quantify this, 
we use the Delaunay triangulation to determine the near- 
est neighbor network of our particle configurations. The 
nearest neighbor network tells us the number of near- 
est neighbors, or coordination number, of each particle. 
For a system of particles in a periodic plane, the aver- 
age coordination number is always si»i£. Particles in a 
triangularly ordered region will be six-coordinated, while 
disruptions in the lattice will lead to particles with coor- 
dination numbers greater than or less than six. A defect 
is defined as any coordination number other than six. 
These non-six-coordinated atoms may be thought of as 
disclinations of charge n, their coordination number be- 
ing 6 + n. 

The most common type of disruption, or defect, is a 
five- or seven-coordinated particle. These may be inter- 
preted as disclinations of charge plus or minus one. Two 
oppositely charged disclinations may be thought of as a 
dislocation. More complex arrangements of disclinations 
are possible, such as dislocation pairs and grain bound- 
ary loops, but in our analysis we have only considered 
individual defects. The defect fraction, fd = 1 — Nq/N, 
is defined as the fraction of particles that do not have six 
neighbors, where N is the number of particles in the sys- 
tem, and Nq is the number of six-coordinated particles in 
the system. Remembering that dislocations are made of 
two bound disclinations of opposite charge, and that dis- 
locations become unbound above the melting point, we 
can expect the defect fraction to experience a jump at 
the melting points. Additionally, at low temperatures 
we can expect an energy gap to occur, which is the en- 
ergy cost to create a dislocation pair. Because the overall 
disclinicity of the system must be zero, as well as the net 
Burgers vector of any dislocations, the lowest-energy de- 
fect excitation is a dislocation pair of opposite Burgers 
vectors. In practice this is usually two pairs of 5- and 
7-coordinated particles. This leads to an exponential be- 
havior in the defect fraction, fd — e~^ A , where A is the 
lowest-energy for a defect type excitation of the system. 




5 10 15 20 25 30 35 40 45 



FIG. 1: The Delaunay triangulation for N = 1600 particles 
at T=0.70. Defects are shown in red. 
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FIG. 2: The Delaunay triangulation for N = 1600 particles 
at T=0.90. Defects are shown in red. 



B. Unbinding of defects 

In Figures I1I2I3I the Delaunay triangulated configura- 
tion of a 1600 particle system is shown at temperatures 
0.7, 0.9 and 1.1 respectively. The defects are shown in 
red. At low temperature as demonstrated in Figure [TJ we 
see that defects occur in quadruplets consisting of two 5- 
coordinated and two 7-coordinated particles. As the tem- 
perature is raised to 0.9 (Figure [5]) we can see isolated 
dislocations (one 5-fold coordinated atom bound to a 7- 
fold coordinated atom). At yet higher temperature, such 
as 1.1 (Figure [3]) we can observe isolated disclinations. 

This can also be seen in the pair distribution func- 
tions 377 (r), 5(55 (r), and #57 (r), for pairs of 7-coordinated 
particles, pairs of 5-fold coordinated atoms and for 5- 
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FIG. 4: The pair distribution function for 7-coordinated par- 
ticles, 377(7"). The peak for T=0.70 extends to ~ 50. 



fold- 7-fold coordinated atoms respectively. In Figure 21 
a sharp peak in 377(7") is observed at low tempera- 
tures (T = 0.70), indicating that dislocations are tightly 
bound. At higher temperatures (T = 0.90 and T — 1.10), 
the peak in 377(7-) is greatly diminished, and dislocations 
become first weakly bound (T = 0.90) and then com- 
pletely unbound (T = 1.10). 955(f), while not shown, 
behaves qualitatively similar to .977(7"), as both are rep- 
resentative of the pair distribution of dislocations. 

The pair distribution function for disclinations, 357(7"), 
is shown in Figure [3] While the sharp peak at low 
(T = 0.70) and intermediate (T = 0.90) temperature is 



expected, the peak at T — 1.10, while quite lower, is still 
very substantial. This indicates that disclinations have 
not become completely unbound, and indeed it is difficult 
to find isolated disclinations in the snapshot configura- 
tions presented in Figure [3] When isolated disclinations 
do occur, they are still next-nearest neighbors with at 
least one other disclination of opposite charge. 



C. Defect fraction 

According to the KTHNY theory, disclinations remain 
very tightly bound below T m . Above T m , the disclina- 
tions are screened from one another by the presence of 
free dislocations yet remain bound, albeit by a weaker 
logarithmic binding 3 . Thus, we expect a proliferation of 
defects to occur around T m , and to continue growing un- 
til somewhere above Tj, where a saturation should occur. 
In Figure [6] we show the average defect fraction as a func- 
tion of temperature. At low temperature, there are very 
few defects, while at high temperature there is a consid- 
erable fraction of the system that is defected. In between, 
there is a region of rapidly increasing defect fraction, from 
T = 0.8 to T = 1.0. This can be quantitatively verified 
by calculating the temperature derivative of the defect 
fraction, which is indeed found to have a broad peak in 
this temperature region. The overall shape of dfd{T)/dT 
is very similar to that of the specific heat capacity, to be 
shown next. Additionally, we can see some size depen- 
dence in the region 0.6 < T < 1.0, although this seems to 
be an issue mostly for comparisons of the smallest system 
size (N = 1600) to the larger system sizes. 

The specific heat per particle at constant volume, cy, 
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FIG. 6: Fraction of defects, fd, as defined by the fraction of 
non-six-coordinated particles in the Delaunay triangulation, 
fd = 1 — Ne/N. The rapid rise in fd from near zero to almost 
25% is a possible sign that dislocation and/or disclination 
unbinding is occurring. 




FIG. 8: The distribution function shows ordering at low tem- 
peratures, as shown above for T = 0.50, while at higher tem- 
peratures, such as T=2.00 shown above, there is a loss of 
order over moderate length scales. 
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in the specific heat at both T m and T^. However, it is 
not clear whether this will be visible above background 
contributions to the specific heat. Either way, the peak 
in specific heat points to a rearrangement of order in the 
systems studied. Also, if we look at the distribution func- 
tion (Figure [8]), we see ordering at low temperatures, and 
fluid behavior at high temperatures. Overall, it is clear 
that there is a phase transition occurring, with a disor- 
dered fluid state at high temperatures and an ordered 
state at low temperatures. 



FIG. 7: The presence of a peak in the specific heat is indica- 
tive of a phase transition. Interestingly, the peak near T = 0.9 
appears to lessen in magnitude as the system size is increased. 

D. Defect excitation energy 



can be calculated from the energy fluctuations, 
1 (E 2 ) - (E) 2 



c v 



N k B T 2 



(2) 



where E is the total energy of an N particle system. We 
have calculated the specific heat and show it as a func- 
tion of temperature in Figure [7] One can see a broad 
peak in the specific heat per particle. According to the 
KTHNY theory, there should be an essential singularity 



In the KTHNY theory, dislocations are bound at low 
temperatures, and there is a defect core energy associ- 
ated with their creation. This leads to an energy gap, and 
thus using the Arrhenius law, we expect fd = e -2E °/ fcBT , 
where we have used 2E C because dislocation pairs are 
the lowest energy excitation (isolated dislocations are 
forbidden). In Table Q] we show the defect activation 
energy as calculated by the Arrhenius law at low tem- 
peratures. Taking the low temperature limit, we find 
E c = 1.49 ±0.01. 
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IV. ORDER PARAMETERS 

A. Definition and temperature dependence 

Let us define a global order parameter of translational 
order, 

, N 

V& <5 = ^E ex p(^-^), (3) 

J'=l 

where G is a reciprocal lattice vector, and fj is the posi- 
tion vector of particle j. If there is translational ordering 
in a system, then clearly will be non-zero if G is 
a reciprocal lattice vector of the appropriate lattice ge- 
ometry. Due to the shape of our simulational cell, at 
low temperatures this will be a triangular lattice with 
nearest-neighbors in the x-direction. At high tempera- 
tures, no translational ordering is present, and all possi- 
ble values of G should give the same (qualitative) result. 
However, at intermediate temperatures, it may be possi- 
ble for there to be some degree of translational ordering 
that is not strictly commensurate with our simulation 
cell. Indeed, we have observed "canted" solid phases at 
intermediate temperatures, where we find partial trian- 
gular order with nearest neighbors in a direction titled 
from the x-axis by a small angle. In this case, if G for 
the triangular order commensurate with our simulation 
cell is used, tyg will be found to be zero. However, if 

we use an appropriate G for the order present, \& g will 
be found to be non-zero. For this reason, we define the 
true translational order to be the maximum value of ^ g 

for all G. In practice, it is not possible to perform this 
optimization for each Monte Carlo configuration, so we 
make the following assumptions. First, due to the nature 
of ordering in two dimensions, we assume any lattice will 
be triangular. Second, because the density of particles 
is fixed, we assume the lattice spacing in said triangular 
solid to be the same as that for the commensurate cell. 
Thus, we keep the magnitude of G constant, and simply 
determine the direction of solid ordering for each config- 
uration by looking at the average bond direction between 
nearest neighbor particles. This turns out to be a good 



Temperature 


N=1600 


N=6400 


N=25600 


0.50 


1.49946(30) 


1.4919(19) 


1.4873(14) 


0.55 


1.50267(85) 


1.4884(31) 


1.4778(22) 


0.60 


1.4996(14) 


1.4635(43) 


1.4252(31) 


0.65 


1.4872(26) 


1.3908(59) 


1.3480(29) 


0.70 


1.4543(30) 


1.2795(49) 


1.2835(15) 


0.75 


1.3640(54) 


1.2027(19) 


1.2390(26) 



TABLE I: Defect activation energy for various system sizes 
and temperatures, as computed using the Arrhenius law. The 
numbers in parentheses are the uncertainty of the trailing 
digits. 



estimate of the true translational order for a system, but 
it must be remembered that it is strictly speaking a lower 
bound. 

The local order parameter which measures the degree 
of six-fold orientational ordering is defined as 

n(i) 

* 6 (n) = 7TE^ (4) 

nit) 
V > j= i 

where Otj is the angle of the bond between particles i 
and j and the sum over j extends over all n(i) near- 
est neighboring atoms found by the Delaunay triangula- 
tion. The global order parameter associated with bond- 
orientational order is obtained as an average over all par- 
ticles. 

1 N 

*e=jvX>(3)- (5) 

In a perfectly bond-ordered triangular solid, we have that 
n(i) = 6 and 6{j = 7r/3 for all j = 1, .., 6. In such case 
1(^6)1 = 1- In the low temperature phase, there is bond- 
orientational order, so (^e) should be a point on the 
perimeter of a circle with a radius approaching unity as 
T — > 0. In the hexatic phase there is quasi-long-range 
bond-orientational order, which implies that the distri- 
bution of (^e) should become a ring in the imaginary 
plane. In the isotropic phase both (^q) and (^q) should 
be distributed around zero value. 

In the top panel of Figure [§] we show the second mo- 
ment of the translational order parameter, "Jg . There 
appears to be a transition from a translationally ordered 
phase at low temperatures to a disordered phase at higher 
temperatures. In the ordered phase there is a clear rela- 
tion between and system size. We will explore this 
relation in a later section, but for now let us point out 
that this finite-size scaling relation begins to break down 
above T = 0.60. This is expected within the KTHNY 
theory of melting due to the unbinding of dislocations. 
However, on closer inspection, the behavior of the curves 
for N = 6400 and N = 25600 in the region 0.6 < T < 0.8 
is not a smooth connection of the data at higher and 
lower temperature. This is due to our measured quantity 
being a lower bound of translational order. 

Also shown in Figure H] is the second moment of the 
bond orientational order parameter, (bottom panel). 
At low temperatures there is substantial bond orienta- 
tional order. Below T = 0.70 there is very little de- 
pendence of on system size. As the temperature is 
increased, ^ begins to show a marked dependence on 
system size as well as a steep decline in value as we ap- 
proach the high temperature disordered phase. 

B. Order parameter distribution 

The main prediction of Halperin and Nelson^ and 
Young 4 is that if two dimensional melting is the result 
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FIG. 9: The second moment of the translational (top) and 
bond orientational (bottom) order parameters. 
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FIG. 10: Intensity plots of the probability of (left column) 
and *l/6 (right column) on the complex plane for T=0.70 (top 
row), T=0.90 (middle row), and T=1.10 (bottom row). 



der. This is exactly what is expected of the hexatic fluid 
phase. Finally, at T = 1.10 (bottom row) we see that 
both order parameters are distributed about the origin, 
indicating an isotropic fluid phase of no order. 



of dislocation unbinding, as proposed by Kosterlitz and 
Thoulessi, then a second unbinding transition (of discli- 
nations) is required to reach an isotropic fluid state. This 
implies the presence of a novel hexatic fluid phase. In 
Figure [10] we show an intensity plot of the distribution 
of 'Fg and of 'Fg on the complex plane for three different 
temperatures. 

At T = 0.70 (top row), our calculation of the distri- 
bution of the order parameters finds a ring of values for 
^q, while 'Fg is localized in a small region away from 
the origin (a very narrow peak showing as a "star" along 
the positive real axis). This is consistent with the pres- 
ence of long-range bond orientational order (K^e)! > 0), 
while the ring of >Fg values is expected for quasi-long- 
range translational order. At T — 0.90 (middle row), we 
that is clustered about the origin, indicating a lack 
of translational order. Interestingly, >Fg now shows a ring 
of values about the origin, indicating quasi-long-range or- 



V. CRITICAL EXPONENTS 

In the topological solid phase, the scaling form for the 
second moment of the translational order parameter is 
~ L~ n , where L is the (linear) system size and r\ is 
a critical exponent. In the hexatic fluid phase, a similar 
relation holds for bond orientational order, (\&g) ~ L~' ,e . 
According to the KTHNY theory of melting, the critical 
exponents r\ and ?yg will have specific values at melting. 
The translational critical exponent is bounded at lower 
melting temperature: 1/4 < rj{T m ) < 1/3. Additionally, 
the bond orientational critical exponent grows from zero 
at T m to 1/4 at Tj, and is related to the translational 
correlation length: r) 6 (T) ~ £(T)~ 2 £. 

By plotting (or (^@)) versus L on a log-log plot, 

we can find r\ (or ?yg). To demonstrate the validity of this 
scaling law and that our results are not limited by system 
size, in Figure [TT] we plot (^g) versus InL for all system 
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FIG. 11: Scaling of < > with system size L, shown on 
a logarithmic plot. Results for T=0.84 are shown as blue 
circles, while data collected at T=0.92 is represented by red 
squares. In both cases, the data for the three smaller systems 
was fit to the equation In < *I/g >= —r]6 In L + const, and 
the result is plotted as the dotted and dashed lines (the solid 
line is the KTHNY value of tjq at Ti). In both cases, the 
value of In < *6 > of the largest system size (N = 102400) is 
reasonably close to the value expected from scaling. 



sizes at the two temperatures where we have results for 
the N = 102400 system. Results of linear least squares 
fits to the three smallest system sizes (used to gener- 
ate the data for Figure IT2|) are shown as a dotted blue 
line (for data at T=0.80), a dashed red line (for data 
at T=0.84), and a long-dashed green line (for data at 
T=0.92). For the higher temperature, the N = 102400 
data falls directly on this line, within error bars. At 
T=0.84, however, the N = 102400 data indicates that a 
smaller value for 776 at this temperature may be necessary. 
This could either be due to the (presumably) large trans- 
lational correlation lengths at this temperature, which 
would invalidate results for small system sizes, or per- 
haps a very long relaxation time. Either way, from our 
data it is clear that by T=0.92 the KTHNY value of i] e 
at Ti has been well passed. 

In Figure [12] we show the extracted values of 77 and 
t?6, the critical exponents of translational and bond ori- 
entational order. In both panels, we show our results 
as red circles. In the top panel, we can see that r\ 
crosses the KTHNY melting value in the temperature 
range 0.6 < T < 0.65. In the bottom panel we show the 
critical exponent of bond orientational order, 776. This 
exponent crosses the KTHNY melting value (see dashed 
line) at a temperature near 0.89, in close agreement with 
the value for Tj derived from the divergence of the corre- 
lation length ^6 obtained in the next section. This value 
for Ti is also in good agreement with the value reported 




FIG. 12: Anomalous dimensionality of (top) the translational 
order parameter and (bottom) the bond orientational order 
parameter. Our current results are shown in both figures as 
red circles. Shown for comparison are the results of Udink 
and van der Elsken— (blue triangles, both figures). In the 
top figure, the dashed and dotted lines represent the lower 
and upper bounds of r\ at T m , according to KTHNY theory; 
in the bottom figure, the dashed line represents the predicted 
value of 7j6 at T . 



in Ref. 20. However, in Figure [12] we also show the alge- 
braic exponents reported by Udink and van der Elsken 15 . 
In both panels, we can see that their values cross the 
KTHNY melting zone at higher temperatures than our 
values. We believe this disagreement may be due to insuf- 
ficient thermalization time in their study, as this could 
lead to artificially low values of the critical exponents. 
We should sound a note of caution here in regards to the 
scaling of (*&q). Because our measurements for "F^, are 
lower bounds, it is possible that the extracted exponents 
fy(T) are not correct in the temperature regime where G 
is no longer commensurate with the simulation cell, as is 
the case for T > 0.6. 



VI. CORRELATION FUNCTIONS 

The correlation function for bond orientational order 
is given by 



C 6 (r) =< MrW 6 (0) >, 



(6) 



where ip&{r) is the local bond-orientational order param- 
eter defined in Sec. [IV] In the isotropic fluid phase, the 
asymptotic form of Cs(r) is ~ exp(r/£e)^. At shorter 



9 



distances, however, a power law decay comes into play, 
such that as £ 6 diverges as 7$ is approached from above, 
then the asymptotic form becomes Cq{t) ~ r~ m at Ti 
and below, with r}%{Ti) — 1/4. Additionally, we observe 
oscillations in Cq (r) that seem to decay with an exponen- 
tial envelope. Thus, we used the following fitting form for 
the bond orientational correlation function for distances 
r much less than the system size L: 



C 6 (r)=A- 



-r/te 



r'hi 



+ B sm(kr + 5)- 



-r/£ 



r'l 



(7) 



We have used a particle-centric definition of the bond 
orientational correlation function, so in our calculations 
of Ce(r) there will be an influence from g(r), the pair 
distribution function. In the limit of perfect bond ori- 
entational ordering (tpe — 1 everywhere), Cq{t) and g(r) 
will be equivalent. We approximate the oscillatory por- 
tion of Cq (r) which is due to the translational atomic ar- 
rangement using a damped oscillator. The periodic form 
is captured by using sin(fcr + 6), where k is expected 
to be near the first reciprocal lattice vector in magni- 
tude (~ 6) and S is just a phase-shift parameter. The 
size of the oscillations are expected to decay exponen- 
tially in the fluid phase, and algebraically in the hexatic 
phase, so we add also a power law, ending up with a 
term sin(kr + 5)r~ v e~ r . An example fit is shown in 
Figure [T5] Note that the fitting procedure returns pa- 
rameters much more precise than the error bars in Fig- 
ure [13] would indicate are possible. This is due to the 
high degree of correlation between neighboring points of 
Ce(r). In fact, up to a separation of 3 the values of Cq{t) 
are still 99% correlated! This simply means that the rela- 
tive form (including the rate of decay) of C§{r) is consis- 
tent between our various calculations, remembering that 
we average the values of 100 independent parallel Monte 
Carlo simulations. 

We wish to note that we observe an upturn in C§(r) 
as r approaches L/2. At temperatures closer to melt- 
ing (larger correlation lengths), the upturn occurs fur- 
ther from L/2. Next, we would like to determine a char- 
acteristic distance R for a given finite-system of linear 
dimension L so as to stay away from this upturn due to 
finite-size effects. Namely, we wish to limit the range of r 
in our fit of the correlation function to the form given by 
Eq. [7] in the range & < r < R. Let us assume a periodic 
form for the correlation function: 



C (r) = A ( ex P(~ r l&) , exp{-(L-r)/^ e ) 



(8) 



Neglecting the power-law term, the upturn is expected 
to occur when the L — r terms are a significant fraction 
of the r terms. Thus, 



(9) 



where R is the distance at which the L — r terms are a 
fraction x of the r terms. Using x = 0.05 or 5%, this 




FIG. 13: Example of fitting the bond orientational correlation 
function, C@, to the form shown in Equation [T) The data is 
for N = 25600 particles at T = 0.97. The critical exponents 
are fixed at their maximum values, r) = 0.33 and rjs = 0.25. 
The extracted correlation lengths are £ = 7.40 ± 0.19 and 
£e = 32.6 ±0.7. 




FIG. 14: Correlation lengths of the bond orientational order 
parameter as determined by fitting the bond orientational cor- 
relation function to the form mentioned in the text. 



leads to 



2 2 



(10) 



In Figure Q3] we show ^6 {T) as determined by fitting 
Gs(r) in the range < r < R. These values were fit to 
the KTHNY form of the expected divergence of ^6 as 
is approached from above: £e(T) — Aexp(6/t I/ ), where 
t = (T - Ti)/Ti and v = 1/2. This fit gives a value for 
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FIG. 15: Demonstration of finite-size-scaling by plotting the 
dimensionless quantity Z/* 76 versus ln(L/^6), i.e., the log- 
arithm of the ratio of the finite-system-size to the correlation 
length, for various size systems. 




FIG. 16: Example of fitting the pair distribution function, 
g(r), to the form shown in Equation 1131 The data is for 
N = 25600 particles at T = 0.97. The critical exponent r\ is 
fixed at its maximum value, 1/3. The extracted correlation 
length is £ = 8.09 ± 0.04. 



Ti near 0.89. 

Using the calculated correlation length and critical ex- 
ponent in Figure [TS] we plot the dimensionless quan- 
tity L t,6 (^S>q} as a function of the dimensionless ratio 
ln(£/^e) for all size-lattice considered here. Notice that 
the data collapse onto the same scaling function using 
the same values of the parameters for our fit to £e(T) 
shown in Figure [HI and also setting r] 6 — = 1/4. 

This provides additional support for the theory. 

VII. DISTRIBUTION FUNCTIONS 

In the disordered phase (T > Tj) the distribution func- 
tion can be obtained as an angular average of the bond- 
orientational correlation function Cg(r) as 

g(r) = l + ±- e i6 - p C 3 {r)d<t>. (11) 
27r Jo 

The integration of e r will give us a zeroth-order Bessel 
function of the first kind, Jq(Gt), and using the KTHNY 
form of the translational correlation function, Cg[f) ~ 
cxp(— r/^)r - '', we wind up with the following form for 
the radial pair distribution function (in the high temper- 
ature limit): 

g(r^oo) = l + AJ Q {Gr)e- r l ( >r- r ' (12) 

where A is some amplitude. 

The G that we use here is the same as in the definition 
of the translational order parameter, namely we use the 



first reciprocal lattice vector of the idealized triangular 
lattice that is commensurate with our simulation cell. For 
the density considered (pa 2 — 0.873), this means G ~ 
6.3c~ 2 . Thus Gr is quite large for moderate values of r, 
and we can use the asymptotic expansion of Jo, namely 
Jo(x — !• oo) = a/2/ttx cos(x — 7r/4). Thus in practice we 
fit g(r) in the disordered phase to the following form, 

e -r/€ 

g (r^^) = l + Acos(kr + S)-^ TT ^. (13) 

An example fit is shown in Figure 1161 In Figure [T7] we 
show the correlation length of translational order as cal- 
culated by fitting g(r) to the above form. Results are 
shown for the N = 25600 and N = 102400 particle sys- 
tems. Clearly, £ remains finite even as the orientational 
correlation length diverges. However, there are some dis- 
crepancies in our values of £. At T — 0.92, the value of £ 
extracted from the N = 25600 particle system does not 
agree with the value for N — 102400 particles. 

We believe that some of this difference may be at- 
tributable to finite size effects. Additionally, there is also 
the possibility that the N = 102400 particle system has 
not fully thermalized. While we have tried to ensure that 
the data for this largest system is completely thermalized, 
it can be very difficult to distinguish between stable and 
metastable states. In either case, we can not consistently 
fit all the data to the KTHNY form, £ = Aexp(b/f ), 
so instead we have made the fit for only the N = 25600 
data. The result of a fit with T m = 0.61, A = 0.00311 
and B = 6.62 using v — 0.36963 is shown in Figure [T71 
as the red curve. In addition, a few other curves are also 
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FIG. 17: Correlation lengths of the translational order pa- 
rameter as determined by fitting the pair distribution func- 
tion to the form mentioned in the text. The range of the fit 
is from 2£ to 4£, with r\ fixed at its maximum value of 1/3. 
At T = 0.80, the fitting range is from £ to 2£. The solid lines 
are obtained from the KTHNY form, £ = A exp(6/t 1/ ), using 
v = 0.36963 and various values of the other parameters. Our 
best fit corresponds to the red curve. 



shown for different values of these parameters with the 
same value of T m = 0.61 whose significance is discussed 
next. 

In Figure [18] we show the approximate validity of 
finite-size scaling by plotting the dimensionless quan- 
tity L^^q) versus ln(L/£), i.e., the logarithm of the 
ratio of the finite-system-size to the measured correla- 
tion length, for two different size systems using the lower 
bound of j] = 0.25 according to the KTHNY theory, 
namely 1/4 < 77 < 1/3. The best collapse is obtained for 
the parameters A = 0.02192, B = 4.89 and T m = 0.61 
shown as the top curve in Figure [T5] (Note that we have 
used the constraint T m > 0.6 as indicated by the be- 
havior of the critical exponent 77). Using the values of 
the parameters obtained for this "best" collapse we ob- 
tain the curve for £(T) shown as a green line in Fig. 1171 
The collapse obtained by using the parameters obtained 
by the best fit to the correlation length (corresponding to 
the red curve in Fig. IT7l) is shown as the graph at the bot- 
tom. We have also included two more fits of both types 
of data, obtained by using parameter values between the 
above two extremes. We can observe that while we do 
not obtain the best fit of both sets of data (i.e., collapse 
of L ri (^ 2 G ) versus ln(L/£) (Tig. [T5|) . and the temperature 
dependence of £ (T) (Fig. [TTJ) for the same values of these 
parameters we see that the values of T m and b are close, 
only the prefactor A cannot be accurately determined. 
We feel that the overall quality of fit is reasonable given 
the fact that we had the difficulty discussed above in de- 
termining the correlation length associated with transla- 
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FIG. 18: Demonstration of finite-size-scaling by plotting the 
dimensionless quantity L^^Sq) versus ln(L/£), i.e., the log- 
arithm of the ratio of the finite-system-size to the measured 
correlation length, for the two size systems, using r\ — 1/4 and 
for various parameters. Due to our calculation of (9q) being 
a lower bound translational order, only data for T > 0.8 are 
shown. 



tional order. 

Several experimental investigations 21,22 have used the 
decay of the envelope of g{r) to extract £. The resulting 
values of £ appear not to diverge across the melting tran- 
sition, so perhaps there is some shortfall in using g(r) 
to get £ at low temperature. For instance, Murray and 
Van Winkle observe a finite peak in £, while for ^6 a 
divergence is seen to occur—. 

Regardless of these differences, if we plot the results 
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T 

FIG. 20: Simplified 20 Binder ratio of the bond orientational 
FIG. 19: Correlation lengths of the translational order pa- order parameter (lines are guides for the eye). The inset shows 
rameter are compared to £e(T). U(^e) near the crossing temperature, T c . 



for £(T) on the same plot with the results for ^(T) as 
shown in Figure Q1J1 we see clearly that these two correla- 
tion length are very different and the differences between 
these various fitting forms for £ are not significant on this 
scale. 



VIII. BINDER RATIOS 

A central concept in finite-size scaling theory is that 
any dimensionless quantity should be a function of di- 
mensionless ratios of the finite-size length (L) of the sys- 
tem to the correlation length £(T) which emerges natu- 
rally and it diverges near the critical point 23 . Therefore, 
close enough to the critical point a dimensionless quan- 
tity becomes a scaling function f(L/£). At precisely the 
critical point, where the correlation length diverges, all 
dimensionless quantities are expected to be independent 
of the system size. 

A straightforward way to construct a dimensionless 
variable is to take the ratio of cumulants. A simple non- 
trivial ratio is the so-called Binder ratio^ of the fourth 
and second cumulants, 



U(x) = 1 - 



((*-(*)) 4 ) 

2\2 ' 



3<(a 



(14) 



As mentioned above, this( dimensionless) variable is ex- 
pected to be system-size independent at a critical point. 
Hence, if the values of U{x) for several system sizes are 
plotted across a continuous phase transition, they should 
cross at the critical point. This is the standard way of 
estimating for example the critical temperature of a ther- 
mal phase transition using the method of Binder ratios. 



In the case of melting in two dimensions, however, we 
have seen that there are two correlation lengths: one for 
translational order, and another for bond orientational 
order. Clearly, if we approach very close to either Tj or Tj 
only one of these two correlation lengths dominates. For 
example if we approach T m sufficiently close from above, 
£ becomes very large and £g is infinite. Thus, there is 
only one finite correlation length. When we approach 
Ti from above, both £ 6 and £ are finite, but if we are 
sufficiently close to Tj, & 3> £ and so we can neglect the 
influence of £. In practice, however, because grows 
very rapidly as the temperature Tj is approached and we 
can only study finite-size size systems, the size of £ is not 
necessarily negligible as compared to the size of . This 
implies that the scaling function becomes /(L/£, L/^q). 
As we have shown in the previous section, £ is still finite 
when £g diverges at the upper critical temperature, Tj. 
Thus, the Binder ratio would only be expected to have 
a crossing at Tj if £ <C L, which is not the case for the 
system sizes we have considered (£(T) w 20, half the 
length of the smallest system size). However, depending 
on the exact form of the scaling function, there may still 
be a crossing in the vicinity of Tj. 

Looking at the Binder ratio in Figure [201 we can see 
that there is an apparent crossing of U^q) at T c — 0.92± 
0.01. Although our statistical uncertainty is too great to 
identify the system-size dependence of the Binder ratio 
crossing (see inset), the finite-size scaling theory outlined 
above indicates that T c should approach Tj when £ -C L. 
Thus, while we could use the value of T c as an estimate 
of T , the method obtained earlier for critical exponents 
is expected to yield more accurate results to the system 
sizes considered here. 

While we have also calculated E/(\&g), the shortfalls of 
our estimator for translational order in the temperature 
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region 0.6 < T < 0.8 lead to an inconclusive analysis of 
the Binder crossing. 

Lastly, let us point out that the finite-size scaling the- 
ory discussed in this section should be applicable to any 
dimensionless parameter. In Sec. I VII we demonstrated 
the scaling collapse of the quantity L ?,6 (4'|) when plot- 
ted as a function of L/£q. In light of the analysis above, 
it is clear that we have neglected the £ dependence of 
this dimensionless quantity. In Figure [17] we can see that 
for T > 0.95 £ is more or less constant (£ rj 8). But 
as Ti is approached, £ increases more rapidly, such that 
£(Ti) fa 20. This could explain the scatter seen in the 
scaling collapse of L V6 (^g) shown in Figure [TBI 

IX. CONCLUSIONS 

Wc have shown that several key predictions from the 
KTHNY theory of two-stage continuous melting are seen 
in the classical system of Lcnnard-Jones (LJ) particles in 
two dimensions. 

First, using Delaunay triangulation we can define 
disclinations and dislocations and this allows us to in- 
vestigate the role of defects in the 2D melting process. 
We can clearly observe at low temperature that disclina- 
tions of 5-fold coordinates atoms and disclinations of 7- 
fold coordinated atoms are bound into dislocations which 
themselves are bound into dislocation pairs. Near T m we 
begin to see unbound dislocations and at a higher tem- 
perature we begin to observe unbinding of disclinations. 
The derivative with respect to temperature of the total 
defect fraction exhibits a broad peak near T ~ 0.9 very 
similar to the specific heat peak. Near this temperature 
we find that the short-range peak (main peak) of the 
pair distribution function of the 5-fold coordinated atoms 
and that of the 7-fold coordinated atoms greatly dimin- 
ishes. The pair distribution function of 5-fold- 7-fold co- 
ordinated particles also decreases greatly at roughly the 
same temperature. 

We calculated the distribution of the order parameters 
and ^6 on the complex plane. Below T m , we see 
the characteristic "Mexican hat" -like circularly symmet- 



ric distribution for ^>q, i.e., while the magnitude of ^g is 
finite below T m , its phase fluctuates causing the system 
to lose its translational order. The orientational order 
parameter, ^g, however, remains frozen to a particular 
direction below T m because the system is large enough 
to allow, for all practical purposes, for such a sponta- 
neous symmetry breaking. In the temperature range 
T m <T <Ti, the "Mexican-hat" -like distribution of * G 
collapses to a distribution around zero, while the distri- 
bution of ^6 becomes "Mexican- hat" -like. This could 
serve as a textbook description of the hexatic order. For 
T > Ti the distribution of both ^>g and VPe are centered 
around zero value. 

We also calculated the temperature dependence of the 
second moment of the above two order parameters for 
various size systems and from the size-dependence of the 
results we have extracted the anomalous dimensions, i.e., 
the critical exponents rj and r/Q. 

Furthermore, we calculated the correlation functions 
Coir) and C§(r) of the order parameter ^g and \I>6 re- 
spectively. We find that both are controlled by two char- 
acteristic correlation lengths, one is £(T), which charac- 
terizes the decay of the correlation of the atomic positions 
and the other is £e(T), which provides the decay of the 
bond-orientation correlations. We demonstrate that we 
can accurately extract both £ (T) and & (T) . 

We find that the two correlation lengths £e{T) and 
£(T) have very different temperature dependence, each 
diverging as we lower the temperature at two differ- 
ent characteristic critical temperatures Ti and T m re- 
spectively, obtained by fitting the calculated correlation 
length to the forms suggested by KTHNY theory. Fur- 
thermore, using the calculated correlation length and 
critical exponent t]q, we find that the dimensionless quan- 
tity L m (\&g) as a function of the dimensionless ratio 
ln(L/£g) for all size-lattice considered here collapse onto 
the same scaling function. A similar conclusion is also 
reached for the finite-size scaling of the corresponding 
quantities related to the translational order, i.e., U 1 ^ 2 -) 
versus ln(L/£). This provides additional support for the 
KTHNY theory. 
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